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Abstract. We present a detailed study of the large-scale collective prop¬ 
erties of self-propelled particles (SPPs) moving in two-dimensional het¬ 
erogeneous space. The impact of spatial heterogeneities on the ordered, 
collectively moving phase is investigated. We show that for strong 
enough spatial heterogeneity, the well-documented high-density, high- 
ordered propagating bands that emerge in homogeneous space disap¬ 
pear. Moreover, the ordered phase does not exhibit long-range order, 
as occurs in homogeneous systems, but rather quasi-long range order: 
i.e. the SPP system becomes disordered in the thermodynamical limit. 
For finite size systems, we Hnd that there is an optimal noise value that 
maximizes order. Interestingly, the system becomes disordered in two 
limits, for high noise values as well as for vanishing noise. This remark¬ 
able hnding strongly suggests the existence of two critical points, in¬ 
stead of only one, associated to the collective motion transition. Density 
fluctuations are consistent with these observations, being higher and 
anomalously strong at the optimal noise, and decreasing and crossing 
over to normal for high and low noise values. Collective properties are 
investigated in static as well as dynamic heterogeneous environments, 
and by changing the symmetry of the velocity alignment mechanism of 
the SPPs. 


1 Introduction 

We understand by an active particle, a particle that is able to convert energy into 
work to self-propel in a dissipative medium. There are several strategies to achieve self¬ 
propulsion. (i) Particles may possess an energy depot and be equipped with a motorQ, 
a scenario representative of many biological systems such as moving bacteria [213] . 
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^ Motility assays represent a particular case, where motors are not attached to the moving 
particles, while the energy depot could be considered as extended over the space [T] 
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insects US], and animals m- (ii) Particles may be able to rectify an external driving 
in order to achieve self-propulsion in a given direction, a situation commonly observed 
in non-living, artificial active particles such as vibration-induced self-propelled rods 
and discs mm, light-induced thermophoretic active particles chem¬ 

ically driven particles |15I16I17I18I19I2T1] . and rollers driven by the Quicke rotation 
effect m- Independently of the strategy exploited by the particles to self-propel, these 
systems are intrinsically out of (thermodynamic) equilibrium 0 and even in absence of 
particle-particle interactions exhibit a non-trivial behavior. For instance, fluctuations 
in the self-propelling mechanism can lead to complex transients in the mean-square 
displacement of the particles [22] as well as anomalous velocity distributions [23]. 


Interestingly, most, if not all, active particle systems found in nature take place, at 
all scales, in the heterogeneous media: from bacterial motion in natural habitats [53], 
such as the gastrointestinal tract and the soil, among other complex environments, to 
the migration of herd of mammals across forests and steppes [7] . Despite this evident 
fact, active matter research has focused almost exclusively, at the experimental and 
theoretical level, on homogeneous active systems [25l2fil27l23l28j . Non-equilibrium, 
large-scale properties of active systems such as long-range order in two-dimensions 
as Vicsek et al. |29j reported in their pioneering paper, the emergence of high-order, 
high-density traveling bands [2011], and the presence of giant number fluctuations 
in ordered phases [26132133] are all non-equilibrium features either predicted or dis¬ 
covered in perfectly homogeneous systems. Here we show that most of these non¬ 
equilibrium features are strongly affected by the presence of spatial heterogeneities. 
Moreover, we show that these properties vanish in strongly heterogeneous media. 
More specifically, we extend previous results [23] on the large-scale collective prop¬ 
erties of interacting self-propelled particles (SPPs) moving at constant speed in an 
heterogeneous space. We model the spatial heterogeneity as a random distribution 
of undesirable areas or “obstacles” that the SPPs avoid. The degree of heterogeneity 
is controlled by the average density po of obstacles. We provide numerical evidence 
that indicates that at low densities of heterogeneities po, the SPPs exhibit, below 
a critical noise intensity ryci, long-range order (LRO). For noise intensities p close 
to Pel, the SPPs self-organize, as in homogeneous space, into high-density traveling 
structures called “bands”. We find that as po is increased, bands become less pro¬ 
nounced to the point that for large enough values of po they are no longer observed. 
Our results indicate that in strongly heterogeneous media, i.e. large values of po, the 
large-scale properties of the system are remarkably different from what we know of 
the Vicsek model [29] in two-dimensional homogeneous media. For instance, orien¬ 
tational order for p < pd is no longer LRO, but rather quasi-long range (QLRO), 
with the system exhibiting the maximum degree of order at an intermediate noise 
value Pm, such that 0 < pM < Pci- Moreover, we provide solid evidence that the 
system becomes disordered as p ^ 0. The numerical data suggests the existence of 
a second critical point pc 2 , with 0 < pc 2 < Pci below which the system is genuinely 
disordered. The disordered phase at low p values is characterized by the presence of 
large, dense moving clusters. We show that the particle number statistics is consis¬ 
tent with these observations: giant number fluctuations (GNF) are high0 near pM 
and decrease as p approaches pc 2 , with GNF becoming weaker as po is increased to 
the point that fluctuations become normal. Finally, we investigate and compare static 
and dynamical heterogeneous media, as well as SPPs with ferromagnetic and nematic 
velocity alignment. We show that in all cases there exists an optimal noise intensity 


^ This is due to the energy consumption involved by the self-propelling mechanism and 
the energy dissipated to the medium 
® The associated GNF exponent adopts its maximum value near pM- 
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Fig. 1. Interaction of a self-propelled particle (blue) with an obstacle (green). Angles are 
given with respect to horizontal axis that is directed from left to right. The moving direction 
given by the angle 6{t) evolves in time from an initial direction at ti to a Hnal direction at 
tf final. The temporal evolution of the particle position, with respect to the obstacles, is 
given by the angle a(t). The initial value of a at ti and its final value at t/ are marked with 
dashed lines. Notice, that after a collision, 9{tf) « ce{tf), the difference A — 9{tf) — a{tf) 
is such that Z\ << 1 for either large values of 7o or large values of Ro- In the hgure, 7o = 1 
and Ro = 1. 


that maximizes the ordering in the systems. The reported results might be of a great 
importance for the design and control of active particles systems. 

The paper is organized as follows. In Section [5] we introduce a general and simple 
model for SPPs in heterogeneous media. The collective motion phase, the associated 
kinetic phase transition and ordering properties of the model are studied in Sec. [31 
Anomalous density fluctuations are addressed in Sec. [H In section O we investigate 
the collective properties of SPPs in a dynamical heterogeneous environment, and 
discuss the effect of the velocity alignment symmetry. We summarize our results in 
Sec. |6l 


2 Formulation of the model 

We consider a system of N self-propelled particles, moving with a constant speed vq 
on a two-dimensional heterogeneous space of linear size L with periodic boundary 
conditions. We express the equations of motion of the i-th particle as: 


±i = VoViOi) 


( 1 ) 


= g{^i) 


lb 




sin [q{0j - Oi)] 


h{yii,9i) + 


( 2 ) 


where the dot denotes temporal derivative, corresponds to the position of the ith 
particle, and 9i is an angle associated to its moving direction. Equation m indi¬ 
cates simply that the particle moves at speed vq in direction V(0i) = (cos^i,sin^i). 
Equation m conveys the dynamics of the moving direction of the particle, which is 
parametrized by the angle 9i. The first term on the right-hand side corresponds to a 
velocity-velocity alignment mechanism acting between neighboring particles as in the 
Vicsek model [29], the second term models the interaction of the f-th particle with 






4 


Kcj incciT’^-orl ^-Vio orli^r^r 



Fig. 2. Simulation snapshots of different phases for po = 2.55 x 10 ® for a system size 
Ni, = 19600 {L — 140). The SPP are represented as black arrows, while the obstacles as green 
dots. The bottom panels correspond to macroscopic phases, while the top panels to zoom 
up regions, where the interaction between the SPPs and the obstacles can be appreciated. 
From left to right: (a) clustered phase, rj — 0.01 and r = 0.58, (b) homogeneous ordered 
phase, rj = 0.3 and r = 0.97, and (c) ordered band phase, p = 0.6 and r = 0.73. 


the spatial heterogeneities, and the third term is an additive noise, where = 0, 

= 6ijS{t — t'), and rj denotes the noise strength. The velocity-velocity 
alignment is characterized by three parameters: its symmetry, given by q and that 
we fix for most of this analysis to be ferromagnetic with <7=1, the interaction radius 
Rb, and the alignment strength jb, where the symbol nb{:x.i) represents the number of 
SPPs at a distance less than or equal to Rb from x^, i.e. the number of neighbors of the 
*-th particle. The function controls the relative weight between alignment (to 

the other particles) and heterogeneity/obstacle avoidance. We tested two possibilities, 
both leading to the same macroscopic behavior: g{xi) = 1 and g{xi) = 1 — 0(no(xi)), 
where 0{x) is Heaviside step function. The latter implies that in the proximity of an 
obstacle, the SP particle focuses on avoiding it, without aligning to the neighbors. 
Results shown here correspond to this definition. Finally, the interaction with the 
spatial heterogeneities, which we refer to as “obstacles” or undesirable areas, is given 
by the function h{xi,9i), defined as: 

h(xi,0i) = —^ sin(afc,i - 6»i), (3) 

if no(xi) > 0, otherwise h(xi) = 0. The term no{xi) represents the number obstacles 
at a distance less than or equal to Ro- The angle ak,i is simply the angle of the polar 
representation of the vector x^ — = Fk^t (cosafey, sina^^i), where yk is the position 

of the fcth-obstacle and Fk^i is the norm of the vector. The position of obstacles, 
with k G [l,7Vo], is random and homogeneously distributed in space. The interaction 
between a SPP and an obstacle is depicted in figure [TJ as the particle approaches the 
obstacle, its trajectory is deflected. While here we focus mainly on “obstacles” whose 
position is fixed in time, we will also address briefly, at the end of the paper, the case 
of moving obstacle, i.e. free to diffuse around the space. 

It is worth analyzing few simple limits of the model given by equations ([T|) and (P|) . 
For 7b = 7o = 0, the equations define a system of non-interacting persistent random 
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Fig. 3. Average values of order parameter r for 10 different initial configuration of obstacles 
for three values of the noise strength rj and po = 0.0325. System size Nh = 10000. Notice that 
the average value of the order parameter does not depend on the particular configuration of 
the randomly placed (static) obstacles. The error bars correspond to the standard deviation 
of the time series of r obtained for each obstacle configuration and set of parameters. 


walkers characterized by a diffusion coefficient = v 1 /rf. With 7o = 0 and 7b > 0 
(or equivalently 7o > 0 and No = 0), the model reduces to a continuous-time ver¬ 
sion [35] of the Vicsek model [29]. For 7o > 0, 7b = 0, and No > 0, the equations 
describe a system of non-interacting active particles moving at constant speed on 
an heterogeneous space, where several interesting non-equilibrium features can be 
observed [36] . At low obstacle density, particles move diffusively with a diffusion co¬ 
efficient that is, interestingly, a non-monotonic function of the obstacle density. It 
reaches a minimum at a given, non-trivial, non-zero obstacle density. On the other 
hand, at high obstacle density and for large enough interaction strength 7 b, sponta¬ 
neous trapping of particles can occur. In this scenario, particles move sub-diffusively 
across the space, spending arbitrary long time in the spontaneously formed traps m- 
Here, we focus on the general scenario where 7 o > 0, 7 b > 0, and No > 0. We reduce 
the parameter space by fixing the following parameters: i?b = i?o = 1, 7b = 7o = !> 
Pb = Nb/L'^ = 1, uo = 1 and a discretization time step to Z\t = 0.1. Notice that our 
main control parameters are the noise intensity 77 and the number of obstacles No or 
equivalently, the obstacles density po = No/L'^- 


3 The order-disorder transition 

The system exhibits three distinct macroscopic phases, for any given obstacle density 
Po > 0, as we move from high to low noise amplitude rj. At high noise values, particles 
are homogeneously distributed in space as a disordered gas of non-interacting parti¬ 
cles. Below a critical noise value 77 ^ 1 , the system undergoes a kinetic phase transition 
from the disordered gas to a (locally) ordered phase. For finite size systems, there 
is a symmetry breaking and the ordered phase implies the existence of a net flux of 
particles in a given direction, or equivalently the existence of a preferred direction 
of motion. For values close to onset of collective motion, i.e. close to rjci, and for 
rather low values of po particles self-organize into high-density, high-order traveling 
structures called “bands”, as illustrated in Fig. He). As 77 is decreased further, getting 
deeper in the ordered phase, bands disappear and we observe an ordered phase where 
particles are roughly homogeneously distributed in space, though anomalously large 
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Fig. 4. (a) Order parameter r vs. noise strength rj for various values of the obstacle density 
Po- (b) Order parameter r vs. obstacle density po for various values of noise strength p. 
System size: Nb = 19600 (L — 140). Notice in (a) that curves for po > 0 exhibit a (local) 
maximum. This implies the existence of optimal noise value that maximizes collective motion. 


density fluctuations are present, Fig. [2Kb). If rj is decreased even further, coming close 
to the noiseless limit, counterintuitively the degree of order drops dramatically and 
particles are organized into densely packed moving clusters, see Fig. UKa), that are 
only weakly correlated due to the constant deflections they experience when running 
into obstacles. 

In order to quantify the degree of order in the system, we use the following order 
parameter: 



where ■ ■)t denotes the average over time0. The definition of r is simply, in complex 
notation, the norm of the average velocity of the particles. Values of r larger than 
zero indicate that there is a net flux of particles in a given direction, while r = 
0 corresponds to a disordered system. It is important to stress that the average 
value of r does not depend on the particular configuration of (static) obstacles, as 
long as the obstacles have been placed randomly in the space. For instance, if we 
compare simulations performed with different obstacle configurations, we obtain the 
same average value of r, as shown in Fig. |31 where the error bars represent the 
measured standard deviation in the corresponding time series of the order parameter 
r. In summary, the value of r, Eq. jT) depends only on the noise amplitude rj and 
obstacle density po- 

Fig. [4] shows the dependency of the order parameter r with respect to the noise 
intensity p for various obstacle densities po, panel (a), and with respect to the obstacle 
density for various noise intensities, panel (b). The curve that corresponds to po = 0 in 
Fig.HKa), black curve, exhibits the behavior that we expect in an homogeneous space, 
i.e. as expected in the classical formulation of the Vicsek model m- This reference 
curve indicates that in an homogeneous system, the maximum order is reached in 
the noiseless limit. Notice that for an homogeneous system, as the noise intensity 
is increased, the order parameter r monotonically decreases until the critical noise 
strength rjci, above which the system is fully disordered and r = 0. Fig. EKa) shows 
that the presence of even a small amount of obstacles leads to a qualitatively different 
picture, with the order parameter r exhibiting a different behavior with rj. All curves 
with Po > 0 reach a maximum at a non-zero value of ry, and all decrease as the noiseless 

Since our initial condition is typically a random distribution of particles over space with 
random moving direction, this average is performed after removing the transient. 
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Fig. 5. Evidence of a first-order phase transition at low obstacle densities, (a) Binder cumu- 
lant G and susceptibility x Si'S function of the noise intensity rj. Notice that around rj ~ 0.8, 
G reaches negative values and x peaks, (b) Histogram n(r) of the order parameter obtained 
from time series of r. (c) Time series of the order parameter r. Notice the flip-flops in the 
time series. System size Ni, = 19600 (L = 140). 


limit is approached. This means that for each value of po, there is an optimal value 
of t] that maximizes collective motion. We refer to this rj value as r]M- On the other 
hand, we observe that if we fix the noise intensity rj and vary the obstacle density 
Po, as displayed in Fig. HDb), the order parameter r monotonically decreases as po is 
increased. Notice that curves corresponding to different noises exhibit a quite distinct 
behavior, for instance compare the curve for t] = 0.3, close to the optimal noise rjM 
for most Po values, with the other curves. 

In the following we divide our quantitative analysis into two statical data sets 
that correspond to low obstacle density and high obstacle density, respectively. We 
show that at quite small obstacle densities, the numerical data is consistent with a 
discontinuous (kinetic) phase transition. The numerical data indicates that as the 
obstacle density increases, the traveling bands become weaker until they disappear. 
We show that once we reach such obstacle densities (i.e. for po > 0.1), the order 
is no longer long-range (LRO) in the ordered phase, but rather quasi-long range 
(QLRO). Our results unambiguously indicate that at such high obstacle densities, as 
we approach the noiseless limit, i.e. when particles self-organize into densely packed 
moving clusters as shown in Fig. [5](a), the system is fully disordered, which suggests 
the existence of a second critical point. 


3.1 Low obstacle density 

To characterize the phase transition to orientational order at low obstacle densities, we 
introduce two additional quantities, the susceptibility x and the Binder cumulant [37] 
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Fig. 6. Finite size analysis at low and high obstacle densities, (a) Scaling of r as func¬ 
tion of the system size Nh for low obstacle densities, here po = 2.55 x 10“®, for vari¬ 
ous values of p (color coded). The solid lines correspond to exponential fittings: r{N(,) — 
Too +C* &cp{—Nb/Nf), with Too, C*, and fitting constants. Notice that for Nt, —>■ oo, the 
order parameter reaches a constant value r{Nh oo) —>■ roo, i.e. long-range order (LRO). 
(b) Scaling of r as function of the system size Nb for high obstacle densities, here p = 0.102, 
for various values of p (color coded). The solid curves correspond to power-law fittings: 
r{Nb) = AN^^, where A and p are fitting constants. We define, see text, quasi-long-range 
order (QLRO) when p < 1/16. On the other hand, p = 1/2 implies that the system is fully 
disordered, (c) Exponent p as function of the noise intensity p for p = 0.102. The diagram 
allows us to define two critical points, indicated by the vertical lines: the vertical line to the 
left corresponds to Pc2, while the other one to pd- In between pc2 < p < Pci the order is 
QLRO. The value 1/16 is indicated by an horizontal red line. 


G, whose definitions are given by: 


X = {r^)t - {r)t , 

(5) 

G - 1 

3{r2)? ’ ■ 

(6) 


We use X to determine precisely the position of the critical point pci and G to estimate 
whether the probability distribution of the order parameter r is unimodal and Gaus¬ 
sian. Notice that the definition of G is directly related to the excess kurtosis. FigjSDa) 
shows both quantities as function of the strength p for one of the smallest obstacle 
densities tested, po = 2.54 x 10“^. The peak of y and the sudden change of behavior 
of G at 77 = 0.8 indicates that the critical point is located at pci = 0.8. The drop of G 
to negative values suggests that the probability distribution of r is bimodal m , as 
confirmed in Fig. [5Kb). This finding is the result of abrupt transitions in the value of 
r along time, often called “flip-flops”, from low values (disordered gas phase) to high 
values (ordered phase) and vice versa, Fig. jSKc). In summary, the statical data close 
to the critical point pci is consistent with a discontinuous (kinetic) phase transition: 
negative values of G, a bimodal distribution, and flip-flops as 77 —> pci- 

Our next step is to determine whether at low obstacle densities the observed order 
for 77 < Pel remains present in the thermodynamical limit. This involves a finite size 
scaling. The goal is to obtain the scaling of the order parameter r with the system 
size. If the system exhibits long-range order (LRO), by increasing the system size, 
while keeping constant the particle density pb and obstacle density po^ we expect r 
to saturate to a non-zero value. As measure of the system size we use Nb B-Fig.lHKa) 
shows that at low obstacles densities, here po = 2.55 x 10 r effectively saturates 


® As measure of system size, instead of Nb, we can use either No or L. Since pb and po are 
constant, knowing either L, Nb, or No, we can determine the other two. 




























9 



Fig. 7. Bands, (a) The particle density profile Pl(x) of the bands along the (band) moving 
direction. The horizontal line indicates the density profile corresponding to an homogeneous 
distribution of particles whose value we denote by ph- (b) Maximum pLmax of the density 
profile p(x)l vs . the density of obstacles. The panel shows the quantity pL max!pH — 1, which 
drops to zero for large densities, which indicates the absence of bands in the system. 


with Nb to a non-zero value for a large range of ?; < rjd values 0. Thus, the numerical 
data at very low densities, up to the system size we could reach, is consistent with 
LRO. This implies that in the thermodynamical limit we expect the system to remain 
ordered, i.e. as Ni, — >■ oo, r —>■ Too ( 77 ) > 0 , where Too ( 77 ) is the asymptotic value of r 
in an infinite system, which is a function of 77 , pt, and po- 

At low obstacle densities, as mentioned above, we observe close to the critical 
point Pel a behavior consistent with a first-order (discontinuous) phase transition. 
This seems to be related [3^ with the emergence of high-order traveling bands 0 , 
as shown in Fig. He). Bands are narrow structures that expand through the whole 
system, elongated in the direction perpendicular to their direction of motion. Several 
density profiles, corresponding to various po values, are displayed in Fig. Ha)- To 
construct these profiles, one needs to project the particle positions onto the moving 
direction of the band and make a histogram. Bands exhibit a sharp front and a smooth 
tail. They move across a background gas of disordered (or weakly ordered) particles 
and accumulate particles in the front as they advance forward, while loosing particles 
in the rear. The presence of obstacles strongly affects bands. As po is increased, 
profiles get smoother as illustrated in Fig. H^)- This can be more quantitatively 
seen in Fig. Hb) that shows how the bands vanish for large value of po- In short, 
as the number of obstacles is increased, bands become weaker, with a profile that 
relaxes towards the background gas. At some point, bands and the background gas 
are undistinguishable and bands vanish. This is evident for po > 0.1, where bands are 
no longer observed. 


3.2 High obstacle density 

At high enough obstacle densities, where bands are no longer observed, the system 
behavior is remarkably different. The finite size study reveals that the system is 
unable to reach LRO at such high po values. We find that the order parameter r 
decays with the system size Ni, as a power-law: r oc where the exponent 

® Unfortunately, we can not be sure that this behavior remains for 77 < 0.05. Smaller noises 
are hard to explore numerically. 

^ For instance, when flip-flops are observed, high values of r coincide with the emergence 
of bands. Along the simulation, bands appear and disappear constantly, as flip-flips do. 
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Fig. 8. Number fluctuation statistics, (a) Scaling of the variance of the particle number a{l) 
vs. the average particle number {n)i for different values of noise strength rj for po = 0.102. 
(b) Scaling of g{1) vs {n)i for different values of the density of obstacles po and fixed noise 
strength p = 0.3. In (a) and (b), the scaling g{1) oc {n)\ and g{1) oc - used as reference 

- are indicated by a solid red and a dashed blue curve, respectively. The dashed lines in 
(a) and (b) represent power-law fitting curves, (c) The dependency of the scaling exponent 
P (see text for more details) on the noise strength p at fixed po = 0.102. For po = 0.102, 
Pel ss 0.45 ± 0.05 and po 2 ~ 0.15 ± 0.05, indicated by two dashed vertical lines in the 
figure. In the range pc 2 < p < Pci, the system exhibits Quasi-Long-Ranged Order (QLRO). 
(d) Dependency of the scaling exponent /3 on the density of obstacles at fixed noise strength 
p = 0.3. The region where the system exhibits ordered phases is indicated by a vertical 
dashed line, where order is first LRO and then QLRO. System size Nb = 40000 {L — 200). 


p,{p) is a function of p, as shown in Fig. EKb) for po = 0.102. For noises close to 
the optimal value pM, we obtain exponent values such that 0 < p < 1/16. Here, by 
analogy with Kosterlitz-Thouless transition [5S] we say that when 0 < p, < 1/16 there 
is quasi-long range order (QLRO), while for p > 1/16 we assume that the system is 
disordered. Fig. c) shows that the behavior of p{p) is such that we can define two 
disordered phases, one at high noise values and the other one at low noises, and the 
ordered phase with QLRO at intermediate noises. This implies the existence of two 
critical points, which obey p{pci) = 1/16, with i = 1,2 such that pc 2 < Pci, see 
Fig. EKc). The system exhibits QLRO when pc 2 < p < Pci, while being disordered 
for p < pc 2 and p > pd- Notice that for both disordered phases, p reaches 1/2, in 
particular we observe that ^ 1/2 in two limits: ry —>■ 0 and p —>■ oo. A scaling 

r oc ' corresponds to a fully disordered system with a random distribution of 
moving directions. Interestingly, the densely packed moving cluster phase at low p 
values, as the one observed in Fig. [2Ka), corresponds, at high obstacle densities (i.e. 
for Po > 1), to a fully disordered phase. The presence of only QLRO, or rather the 
absence of LRO, implies that in the thermodynamical limit we expect r —>■ 0 for all 
p > 0 values, i.e. we expect an infinite system to be disordered. Nevertheless, between 
Pc 2 < p < Pel, we expect the SPPs display large correlations in their moving direction 
also in the thermodynamical limit. 
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4 Anomalous density fluctuations and clustering statistics 

A way to characterize the distribution of SPPs over the space is through the study 
of density fluctuations. Particularly useful information is provided by the so-called 
number fluctuations. The idea is to divide the space over which particles move in 
cells of linear size I and count the number of particles in each cell. Let us call n(xi, 1) 
the number of particles in the cell of linear size I whose center is at position x^. We 
are interested in computing the average of this quantity {n{x.i,l))i and its standard 
deviation a. It can be shown m that {n{x.i,l))i = = {n)i and that 

cr(/) = \l{(n{xi,l) - {n)i)f)i = (n)f , (7) 

where the average (.. .)i is performed over the cells the space has been divided into. 
An exponent /3 = 1/2 is expected for a random distribution or particles, while /3 > 
1/2 corresponds to giant number fluctuations (GNF). It has been predicted that in 
spatially homogeneous systems, the SPPs in the ordered phase - far away from the 
band regime - exhibit GNF . This prediction has been observed in simulations of 
SPP particles, where a value of (3 ^ 0.8 was found |88I4()] . Here, we are interested 
in knowing the impact of an heterogeneous environment on the number fluctuations. 
Figure |S] shows how number fluctuations are affected by changing the noise intensity 
77 for a fixed density of obstacles po ~ 0.102 (see Fig. [Ha),(c)) and by varying the 
obstacle density po while keeping the noise fixed, here p = 0.3 (see Fig. |SKb),(d)). 
The top panels of Fig. [5] correspond to the scaling of a{l) with the average number 
{n)i of particles per cell. In the figure, dashed blue lines correspond to a slope 1/2, 
while red solid lines to a slope 1. Undoubtedly, GNF are also present in heterogeneous 
space. Nevertheless, there are important differences with what we know from SPP in 
homogeneous space. For instance, for a fixed (high enough) obstacle density po, GNF 
are suppressed, or at least decrease, as p pc 2 , i.e. (3 adopts smaller values as p 
approaches pc 2 . Fig. Elc). We recall that in homogeneous space GNF are expected to 
be characterized by the same anomalous exponent as 77 —>■ 0. We also point out that 
in both homogeneous and heterogeneous space, number fluctuations become normal 
for p > Pel, i.e. in high-noise disordered phase. This means that at high obstacle 
densities (i.e. for po > 0.1), GNF are stronger - meaning that f3 adopts its largest 
value - at some point in between pc 2 < Pm < Pci, and this seems to occur close to 
Pm- On the other hand, if we fix the noise intensity 77 , we observe that GNF decay 
(i.e., [3 goes down) as the density of obstacles po is increased, reflecting the global 
tendency of the system to go to disorder when po ^ 00 . We find that for po —>• 0, 
f3 0.8 as expected in an homogeneous media, while as po is increased, the exponent 
{3 exhibits two regimes with po, approaching linearly for high obstacle densities 1/2, 
where fluctuations can be considered normal and the system disordered, see Fig.lS^d). 

Another alternative to study how particles are distributed in space is to look 
at the cluster size distribution. As before, we are interested in understanding how 
the presence of obstacles affects the non-equilibrium clustering statistics of the SPPs 
with respect to what we know from homogeneous media |41I42I43] . By “cluster” 
we understand a group of connected particles, such that the distance between two 
connected particles is smaller or equal to the interaction radius. The size or mass of 
a cluster, which we denote here with the letter “m”, is the number of particles the 
cluster contains. Our quantity of interest is the (weighted) cluster size distribution 
(GSD) P{m). Its definition is given by: 

P{m) = lim P{m,t) = lim ^ (g) 

t—¥oo t—^oo 

where n-mit) refers to the number of clusters of mass m that are present in the 
system at time t. The limit is to indicate that we look at the steady state CSD and 
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(a) (b) 



Fig. 9. Clustering statistics, (a) The cluster size distribution (CSD) P(m) for the fixed 
density of obstacles po = 3.25 x 10“^ for different noises r]. The critical noise values are 
? 7 c 2 ~ 0.05 and rjai ~ 0.75. The CSD for r] = 0.9 corresponds to a disordered phase, while 
the other two CSDs to ordered phases. The dashed and dotted-dashed curves provide the 
scaling oc and oc m~°'^ used as reference, (b) P{m) for fixed noise value rj = 0.6 and 

two different values of the obstacle density po. With po = 2.5 x 10“® the system is in an 
ordered phase and bands are observed, while with po = 0.4 the system is fully disordered 
and bands are not observed. The dashed line corresponds to the scaling oc . System 

size Nb = 19600 {L = 140). 


neglect transitory behaviors. Fig. IHDa) shows how the CSD is changed by varying the 
noise rj for fixed po = 3.25 x 10“^. As a reference, the CSD corresponding to a fully 
disordered phase, i.e. p = 0.9, is shown. We find that in between and pM the CSD 
distribution is roughly power-law, P{m) oc m““, with an exponent w ^ 1.18 that falls 
in the range [0.8,4/3] as expected [33]- As we move to lower noise values, e.g. p = 0.3, 
there is a strong depletion of isolated particles and small clusters and particles tend to 
form larger clusters, way larger than those observed close to pci or in the disordered 
phase (notice the log scale in the figure). Fig. [5][b) displays the CSD at fixed noise 
p = 0.6 and two different values of po- For po = 2.5 x 10“^ we observe bands and the 
CSD is again a power-law with exponent w ^ 1.18, similar to what was reported for 
the Vicsek model m in the band regime where uj ~ 1.3 m- The figure evidences 
that by increasing the density of obstacles po, at fixed noise, the functional form 
of CSD is dramatically affected. In particular, it shows that at very large obstacle 
density the CSD becomes exponential as expected for the disordered phase, with a 
well defined average cluster size 0 and not surprisingly the system exhibits normal 
number fluctuations. 


5 Discussion: static vs dynamic heterogeneities and the symmetry 
of the interactions 

There is no reason to believe that static and dynamic heterogeneous environments 
lead to similar large-scale collective effects. In particular, the conclusions drawn from 
the finite size analysis performed with static obstacles, i.e. with a sort of “quenched” 
noise, may not apply to dynamical heterogeneities. One may argue that dynamical 
heterogeneities may be mapped to an effective noise in a SPP system with homoge¬ 
neous space. In this scenario, the dynamical heterogeneities should not affect qualita¬ 
tively the large-scale properties of the system but only have an impact on the critical 
point. A rigorous analysis would require a finite size study of SPPs in dynamical het¬ 
erogeneous environments, which is a very time-demanding numerical task out of the 

® Power-law CSDs may be such that their first moment diverges, while exponential CSDs 
always have a well defined first moment. 












Will be inserted by the editor 


13 


(a) (b) 



Fig. 10. SPP in a dynamical environment, where obstacle diffuse with a diffusion coefficient 
Do- (a) Order parameter r vs. 77 for various obstacles densities po and constant diffusion 
coefficient Do- (b) r vs. 77 for constant density of obstacles po and various diffusion coefficients 
Do - Notice that there exists an optimal noise value even in a dynamical environment. System 
size Nb = 10000 (L = 100 ). 


scope of the current paper. Less ambitious but not less informative, we can analyze 
the impact of a dynamical heterogeneous medium on the collective properties of SPPs 
in a fixed system size. Let us assume that the obstacles now diffusive over the space 
with a diffusion constant Do- The position of the fc-th obstacle obeys: 


ifk = V‘^Do^k(t) 


(9) 


where (Cfc(O) = 9 and ~^')^kj- Fie. ITUl shows the order parameter 

r as function of the (angular) noise 77 , for various values of po and fixed obstacle 
diffusion coefficient Do — 0.7, panel (a), and fixed obstacle density po — 0.102 and 
various obstacle diffusion coefficient Do, panel (b). We find that even for a dynam¬ 
ical heterogeneous environment, there is an optimal noise that maximizes collective 
motion. As the obstacle density is increased, the level of ordering, i.e. r, decreases 
for all angular noises. Fig. ITTir aL On the other hand, we learn that the faster the 
obstacles diffuse, the weaker is the effect of the obstacles. Fig. [TOT b). Moreover, the 
numerical data suggests that in the limit of Do —00 the system behaves again as 
an homogeneous system with its critical point shifted to smaller noise values Q This 
suggests that in this limit effectively the problem can be mapped to an homogeneous 
system with an effective (angular) noise intensity. 

Finally, we may wonder whether the observed optimal value is due to the par¬ 
ticular symmetry of the velocity alignment between the SPPs that has been used, 
i.e. due to g = 1 in Eq. [2l To address this question, we perform simulations with 
SPPs interacting via a nematic velocity alignment, i.e. q = 2 in Eq. [2J that move in 
an heterogeneous medium with static obstacles. Since the alignment is nematic, we 


replace the order parameter by: S 2 = (^ 


where 


)t represents 


Nb l^i=l ^ 

a temporal average after a short transient. For a totally disordered system, S 2 = 0, 
while 5'2 > 0 implies that the system exhibits, for the tested system size, (global) 
nematic order. Fig. 11 shows that optimal noise rjM exists also for SPPs interacting 
via nematic alignment in an heterogeneous environment, as discussed above for fer¬ 
romagnetic velocity alignment. In this case, the optimal noise rjM maximizes nematic 
ordering, i.e. it favors the emergence of a preferred direction of motion, where 50% 
of the SPPs move roughly parallel to it and the other 50% antiparallel to it. Details 
about SPPs with nematic velocity alignment in homogeneous media can be found 


in |:15I44I45| . 


® Compare Fig. llOf bl and Fig. |4{a), curve for po = 0 to see the shift in the critical point. 
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Fig. 11. SPPs with nematic alignment in a (static) heterogeneous medium. Nematic order 
parameter S 2 vs. noise intensity rj for different obstacle densities po- Notice that also SPPs 
with nematic alignment exhibit an optimal noise that maximizes order in the system. System 
size Nb = 10000 {L = 100). 


6 Conclusions 

We have learned that even small levels of heterogeneity lead to qualitative changes 
in the large-scale properties of SPP systems interacting via a velocity alignment 
mechanism. Some of the new statistical features that emerge due to the presence 
of heterogeneities - as for instance the existence of an optimal noise that maximizes 
order [M] - are present in both statical and dynamical heterogeneous media, as well 
as by changing the symmetry of the velocity alignment mechanism of the SPPs. Other 
findings, as the absence of long-range order for high levels of heterogeneity (i.e. for 
Po > 0.1) apply exclusively to static obstacles. In general, we can conclude that in 
heterogeneous environments the physics of SPP systems is different from what we 
know from homogeneous ones, with the presence of obstacles making more difficult 
for the SPPs to spread information about their moving direction across the system. 
From the two information spreading mechanisms |48I49] - the one involving direct 
particle-particle interaction and responsible of order inside individual clusters, and 
the other one involving cluster-cluster information exchange, often through particle 
exchange among clusters - obstacles affect the second oneEj- In particular, due to the 
obstacle presence, clusters get quickly uncorrelated as result of independent collisions 
with the obstacles. At high obstacle densities or low noise amplitudes, particle ex¬ 
change among clusters becomes insufficient to maintain the moving clusters correlated 
and the level of (global) order decreases. 

The spatial arrangement of particles is also strongly affected by the presence 
of obstacles. The high-order traveling bands - reported to emerge in the classical, 
homogeneous, Vicsek model |30l,11] - become less pronounced and even disappear 
as the level of heterogeneities, i.e. obstacles, is increased. At the point where bands 

Information spreading in the context of active particles (without alignment) was studied 
with self-propelled disks that exchange their internal state upon collision [50]. As for align¬ 
ing active particles, there are also two information mechanisms, with particle-particle infor¬ 
mation exchange being way faster than cluster-cluster information exchange. Interestingly, 
fluctuations play a major role even in the absence of cluster-cluster information exchange in 
two dimensions where the physics of the problem is always non-trivial m- 
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are no longer observed, the ordering properties change from long-range to quasi-long 
range, which suggests that in the limit of an infinite system (keeping constant both, 
obstacle and particle density) and for static obstacles, the SPPs cannot maintain 
a coherent migratory route along the (infinite) heterogeneous space: i.e. the system 
becomes disordered. On the other hand, our finite size analysis revealed that at high 
obstacle densities (i.e. po > 0.1), the system exhibits two critical points, one at low and 
another one at high noise value. Finally, the study of density fluctuations indicated 
that the giant-number-fluctuation exponent /3 moves towards 1/2, which corresponds 
to normal density fluctuations, as we approach the noiseless limit, as well as for large 
enough obstacle densities. 

Finally, it is worth mentioning that few experiments with active particle in het¬ 
erogeneous media have been already performed, so far, with active particles without 
alignment: self-propelled janus particles moving on patterned surfaces [46] and speckle 
light fields [47152] . Interestingly, patterned regular environments have been initially 
used to rectify the motion of active swimmers such as bacteria in diluted suspensions 
(i.e. in a non-interacting context) [53154155] . In these systems, volume exclusion ef¬ 
fects and the size of the moving active particles play a central role. Such observations 
have triggered a good deal of theoretical work. For instance, in simulations with self- 
propelled disks (SPD) it has been shown that SPDs can get locally jammed by the 
obstacles [sgB while in simulations with SP rods it has been found that V-shaped 
obstacles can be used to trap particles m- On the other hand, it has been shown 
in simulations with circularly moving active particles that a regular configuration of 
obstacles can be used to filter the active particles [58159] , while narrow channels can 
direct particle motion jSD]. It remains to be seen how the large-scale collective prop¬ 
erties reported here are affected by introducing a velocity alignment mechanism in 
the above mentioned more realistic models. 
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